Below is an example of 20 hypothetical 1ha restoration plots (10,000m2, 100m*100m) distributed at random at three reefs in the Townsville region (Davies Reef, Big Broadhurst Reef, Little Broadhurst Reef). The code takes the Allen Coral atlas geomorphology data (see here for more details), filters all reef polygons larger than 1ha in size, randomly selects 20 polygons across the three reefs, and constructs a 100*100m plot around the centroid of each polygon:
# data packages
library(tidyverse)
library(units)
library(foreach)
# spatial packages
library(sf)
# mapping packages
library(tmap)
library(leaflet)
# create function to get make polygons for plots
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {
# Calculate the coordinates of the rectangular polygon
x <- sf::st_coordinates(input)[1, 1]
y <- sf::st_coordinates(input)[1, 2]
# set parameters
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
sf::st_sfc(crs = 20353)
return(polygon)
}
habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")
townsville_geomorphic <- st_read("/Users/rof011/coraldynamics/data/Davies-Broadhurst-20230906193843/Geomorphic-Map/geomorphic.geojson") %>%
sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
dplyr::group_by(class) %>%
dplyr::mutate(habitat_id = paste0(
gsub(" ", "_", class), "_",
sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
sf::st_make_valid() %>%
# mutate(habitat_id=as.factor(habitat_id)) %>%
dplyr::group_by(habitat_id, class) %>%
summarise(geometry = st_union(geometry)) %>%
mutate(area = round(st_area(geometry))) %>%
filter(class %in% habitats) %>%
drop_na(class) ## Reading layer `Davies Broadhurst' from data source `/Users/rof011/coraldynamics/data/Davies-Broadhurst-20230906193843/Geomorphic-Map/geomorphic.geojson' using driver `GeoJSON'
## Simple feature collection with 2291 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 147.6256 ymin: -18.98455 xmax: 147.8017 ymax: -18.80138
## Geodetic CRS: WGS 84# find restorable habitats
townsville_sheltered_sites <- filter(townsville_geomorphic, class %in% c("Sheltered Reef Slope", "Outer Reef Flat"))
# identify neighbours
neighbors <- st_touches(townsville_geomorphic)
# match pairs
border_rows <- purrr::map2(neighbors, townsville_geomorphic$class, ~ {
any(townsville_geomorphic$class[.x] == "Sheltered Reef Slope") & (.y == "Outer Reef Flat")
})
# Filter rows
indices <- which(unlist(border_rows))
result <- townsville_geomorphic[indices, ]
# combine outer and sheltered
townsville_restorable <- rbind(
filter(townsville_geomorphic, class %in% c("Sheltered Reef Slope")),
result
)
townsville_plots_centroids <- townsville_restorable %>%
filter(area > set_units(10000,"m^2")) %>%
filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>%
drop_na(class) %>%
st_centroid() %>%
st_cast("POINT")
library(foreach)
townsville_plot_outputs <- foreach(i=1:nrow(townsville_plots_centroids), .combine="c") %do% {
tmp <- set_restoration_plot_centres(townsville_plots_centroids$geometry[i], width=100, length=100)
tmp
}Select layers using the layercontrol on the left, use [ ] for full-screen viewing.
library(tmap)
set.seed(1)
townsville_plot_outputs <- townsville_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1)))
townsville_plot_outputs2 <- townsville_plot_outputs[sample(nrow(townsville_plot_outputs), 20), ]
# visualise with thematic maps
map <- tm_view() +
#tm_tiles("Esri.WorldImagery", group = "[Satellite map]", alpha = 0.5) +
#---------- all habitats--------@
tm_shape(townsville_geomorphic |> mutate(class=as.factor(class)), name = "Reef habitats", id="area") +
tm_borders(col = "black", lwd = 0.2) +
tm_fill("class", name = "Benthic habitats", Title = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
"Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
"Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
"Reef Crest" = "coral3"), alpha = 0.6) +
#---------- restorable habitat--------@
tm_shape(townsville_restorable |> mutate(class=as.factor(class)), id="area", name = "Restorable habitats") +
tm_borders(col = "white", lwd = 0) +
tm_fill("class", title="Restorable habitats", name = "Benthic habitats", id="area", palette = c("Sheltered Reef Slope" = "darkslategrey", "Outer Reef Flat" = "darkgoldenrod3"), alpha = 0.6) +
#---------- restoration plots--------@
tm_shape(townsville_plot_outputs2, id="plot_id", name = "Restoration plots (1ha)") +
tm_fill(col="red", alpha=0.6) +
tm_borders(col="black") +
#---------- scale--------@
tm_scale_bar(width=200)
map |> tmap::tmap_leaflet() |>
leaflet::addProviderTiles('Esri.WorldImagery', group = "Satellite map", options=leaflet::providerTileOptions(opacity=0.6, maxNativeZoom=18,maxZoom=100)) |>
# leaflet::addProviderTiles('Esri.WorldTopoMap', group = "<b> [Restoration]</b> habitats", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=10)) |>
leaflet::addLayersControl(position="topleft", overlayGroups=c(
"Satellite map",
"Reef habitats",
"Restorable habitats",
"Restoration plots (1ha)"),
options=leaflet::layersControlOptions(collapsed = FALSE)) |>
leaflet::hideGroup(c("Reef habitats")) |>
leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)